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ABSTRACT 

o , 

■ We model the intermediate time evolution of a "jetted" gamma-ray burst by two 

blobs of matter colliding with the interstellar medium. We follow the hydrodynamical 
^ , evolution of this system numerically and calculate the bremsstrahlung and synchrotron 

images of the remnant. We find that for a burst energy of 10^^ erg the remnant becomes 
spherical after ~ 5000 years when it collects ^ 50Mq of interstellar mass. This result is 
independent of the exact details of the GRB, such as the opening angle. After this time 
^ ■ a gamma-ray burst remnant has an expanding sphere morphology. The similarity to 

a supernova remnant makes it difficult distinguish between the two at this stage. The 
expected number of non-spherical gamma-ray burst remnants is ~ 0.05 per galaxy for 
a beaming factor of 0.01 and a burst energy of 10^^ erg. Our results suggest that that 
I the double-shell object DEM L316 is not a GRB remnant. 
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1. Introduction 



o 

c/3 . If a 7-ray burst (GRB) originate within a galactic disk then the large deposition of energy will 

result in a blast wave whose initial phase produces the afterglow. The late phase of the blast wave 
evolution would result, as noted by Chevalier (1974) in the context of supernova remnants (SNRs), 
^ ' in a cool expanding H i shell. The shell will remain distinct from its surrounding until it has slowed 

[ down to a velocity of ~ lOkms"^. This would happen within 2.3-E'5^^^10^ yr where E^i is the initial 

energy in units of 10^^ erg (Loeb & Perna 1998). The current rate of GRBs is one per ~ 10^ yr 
per galaxy (Schmidt 1999). This leads us to expect a few remnants per galaxy at any given time. 
However SN are 10^ times more frequent and it would be difficult to distinguish between a SNR 
and a GRB remnant. Today there is mounting evidence that some GRBs are beamed (Halpern 
et al. 1999; Harrison et al. 1999; Kulkarni et al. 1999; Sari 1999; Kuulkers et al. 2000; Salmonson 
2000). Beamed GRBs will illuminate only a fraction fj, of the sky, and their rate should be higher 
by a factor of /f^^. With /{, 0.01 we would expect a hundred GRB remnants per galaxy. At an 
early stage the morphology of a "jetted" GRB remnant would be very different from a spherical 
explosion. This can be used to identify those remnants. We study this phase here. 

Both beamed GRBs and SN deposit a comparable (~ 10^^ erg) energy into the ISM. In both 
cases the evolution is expected to be similar since both are in the Sedov (Sedov 1959) regime where 
all the energy is the initial explosion energy and all the mass is in the surrounding ISM. A key 
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distinguishing feature unique to GRB remnants could be their beamed nature which we expect 
would lead to a distinct double shell morphology at intermediate times (the late time behavior of 
the GRB remnant is expected to be spherical in any case). In order to establish how many Hi shells 
are GRB remnants we must find out the expected morphology of GRB remnants and how long they 
remain non-spherical, and distinguishable from SNRs. 

Expanding Hi shells have been found in many spiral galaxies (Tenorio-Tagle &; Bodenheimer 
1988). The interior of these shells is relatively empty and their current expansion velocity is in 
the order of tens of Km s~^. Models for the origin of these Hi supershels involve a large number 
of spatially correlated supernova events (Heiles 1979, 1984), infall of massive gas clouds on to 
the galactic plane (Tenorio-Tagle &; Bodenheimer 1988) and flaring of radio lobes formed by jets 
ejected from the galactic nucleus during an active phase (Gopal-Krishna & Irwin 2000). It has been 
previously noted (Efremov et al. 1998; Loeb &: Perna 1998) that a subset of these Hi supershells 
may be the late signature left by GRBs on the interstellar medium (ISM). Establishing how many 
of the Hi shells are GRB remnants would make it possible to directly estimate the local rate of 
GRBs, determine e, the efficiency of converting the explosion energy into 7-rays, and the beaming 
factor fb (Loeb & Perna 1998). 

We model the intermediate evolution of a beamed GRB by two blobs of dense material moving 
into the ISM in opposite directions and we follow numerically their hydro dynamical evolution. Our 
results can be rescaled to fit a variety of initial energies. We find that at a time of ~ 5000 yr, 
when the ratio Mc^/Eq between the accumulated mass, M, and the initial GRB energy,£'o, reaches 
~ 9.6 X 10^ the remnant becomes spherical, similar in shape to a SNR. This value is independent of 
the exact initial details of our model such as the opening angle, the velocity and the morphology. 
We compare our results to the H I shells DEM L 316 (Williams et al. 1997) previously classified as 
two colliding SNRs. As the accumulated mass there is much larger it is most likely not a GRB 
remnant. 

We describe our model and the numerical methods in section 2 . In section 3 we describe our 
results. Discussion and summary including a comparison with DEM L 316 are given in section 4. 

2. The Simulation 
2.1. Model 

We study the intermediate evolution of the relativistic ejecta that caused a GRB. We assume 
that this matter was ejected in the form of two ultra relativistic blobs moving in opposite directions 
with a bulk Lorenz factor T > 100. The emission from internal collisions in the blobs during an 
early stage comprises the GRB. Late external shocks caused by collisions with circumstellar matter 
produce the afterglow. The matter slows down during this interaction and its bulk Lorenz factor 
r, decreases. The ejecta stays collimated only until F drops below ~ I/Oq, at approximately 
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3.4(E5i,iso/pi)^''^(6'o/0.1)^/^ hr after the GRB, when it has accumulated IQ-^E^i^isoiOo/O.l)'^ M© of 
ISM mass (Sari et al. 1999) where is the initial angular width, £'5i,iso the isotropic energy in units 
of 10^^ erg and pi is the surrounding ISM density in units of lO^^'^ gmcm^"^. Note that in the rest 
of the paper £'51 denotes the actual energy in units of 10^^ erg. At this time the matter will start to 
expand sideways causing, for an adiabatic evolution, an exponential slowing down (Rhoads 1997). 
The ejecta continues to expand sideways at an almost constant radial distance from the source 
Rq ~ 0.3eI(^ p^^^^ pc until it becomes non-relativistic. At this stage, we begin our simulation. 

Without a detailed numerical modeling of the relativistic phase of the ejecta we have only an 
approximate description of how to construct the initial conditions. We expect the angular width of 
the ejecta to be ~ 1 rad. Additionally we are constrained by Newtonian energy conservation which 
yields a relation between Rq, T and Eq: 

i?o - 0.3Ell'p;'/\vo/cr^/' PC . (1) 

The exact shape and energy distribution of the ejecta are uncertain. This may influence the 
numerical coefficient in Eq. (1). However we show that the intermediate and late evolution of the 
ejecta are insensitive to these uncertainties in the initial conditions. This is due to the fact that 
since we are in the Sedov regime, the mass is dominated by the "external" ISM gas which washes 
out any variations in the initial conditions of the ejecta. This same mechanism will smear out any 
non-spherical features in a SN explosion, producing spherical SNRs. 

We realize these initial conditions by two identical blobs moving in opposing directions into 
the ISM. The computational space is a cylinder in which the blobs move along the symmetry axis. 
Both the blobs and the ISM are modeled by a cold 7 = 5/3 ideal gas. The blobs are denser than the 
ISM. We consider various initial densities, angular widths and shapes of the blobs in order to make 
sure that these are indeed unimportant in determining the final morphology. In all runs vq ~ c/3. 
To evolve these initial conditions we use a hydrodynamic code, neglecting radiative effects. 

We note that we have also used the post Newtonian version of our code with the same initial 
blob velocity of c/3. We found that the blobs slow down rapidly and the results were similar to 
the Newtonian results. Within a few timesteps we have a Newtonian system. In essence this is 
equivalent to slightly changing the initial shape of the blobs which, as noted before, is insignificant 
due to the mass dominance of the ISM. 



2.2. The Numerical Method 

The code we use is based on the Newtonian version of the smooth particle hydrodynamics 
(SPH) code introduced in Ayal et al. (2001). The code was adapted for the specific problem at 
hand. The three computationally important features of the problem are the negligible role of 
gravity, the large amount of stationary gas representing the ISM, and the symmetry. The fact that 
gravity is unimportant has made the equations very simple and greatly increased the speed of the 
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code. Implementing single particle time steps allowed us to put in a large volume of stationary gas 
without making a big impact on the computational time. This way most of the computational effort 
is invested in the "important" interacting particles. Typically only 0.1% of the particles actually 
move during each timcstep. 

The initial conditions are set up so that the blobs are on the z axis with a velocity along the z 
axis. With these initial conditions there is a rotational symmetry about the z axis and a reflection 
symmetry about the x — y plane. We implement the reflection symmetry exactly and the rotation 
symmetry only approximately. The implementation, similar to Libersky et al. (1993), consists of 
considering only 1/8 (one quadrant) of the computational volume and adding reflecting boundary 
conditions on the three inner boundaries defined by the x — y, x — z^ and y — z planes. We implement 
no outer boundary condition since the surrounding gas is very cold and there is almost no pressure 
so the outer particles gain only a negligible velocity throughout the course of the simulation. The 
quadrant we evolve consists of 62,800 SPH particles of which only ~ 10 represent the initial blob, 
again stressing the total mass dominance of the ISM. 

The reflecting boundary conditions are implemented using pseudo-particles. At the beginning 
of each timesteps, all particles which intersect one of the boundaries (In SPH each particle has a 
finite size called the smoothing length) are reflected about this boundary and added to the simu- 
lation as additional pseudo-particles. After this is done for all reflecting boundaries, we calculate 
all the time derivatives in the usual manner treating all particles in an equal way. When all time 
derivatives are calculated we delete all the pseudo-particles. Only the "real" particles are then 
advanced in time. This simple algorithm allows us to implement reflecting boundary conditions 
using only a small number (typically 10% of the total number of particles for the three boundaries 
we use) of additional particles. 

The use of SPH requires adding some artificial viscosity in order to resolve shocks. We use 
the standard artificial viscosity (e.g. Monaghan 1992; Benz 1990) consisting of a term analogous to 
bulk viscosity and a Von Ncuman-Richtmycr artificial viscosity term. For the time integration we 
used a second order Runge-Kutta integrator with an adaptive stepsize control. 

3. Results 

We choose the initial velocity to be in the range 0.22c to 0.33c, making relativistic effects small. 
Equation (1) leaves the freedom of choosing two out of the three parameters Eq^ Rq and pi, the 
initial energy, distance and ISM density respectively. In presenting the results we choose £"0 and po- 
To parameterize the evolution of the remnant we utilize the fact that mass scales linearly with the 
initial energy and we define the dimensionless parameter fi = Mc^/Eq where M is the accumulated 
shell mass. We define M as all mass with density above 2po corresponding to the shocked material 
of the shell and to the accumulated ISM mass. All subsequent results are presented as functions 
of 11. Our simulation begins approximately 2 years after the burst, when p ~ 18. In Figure 1, we 
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show the linear scaling of time with //. The scaling relation is t ~ 0.053//(£^5i//Oi)~^/^ yr which is 
close to the Sedov- Taylor (Sedov 1959) blast wave result of t oc R^/'^ oc /x^/^ where R is the radius 
of the blast wave. 

A bow shock forms as each blob collides with the ISM. The shock also propagates in the 
direction perpendicular to the blob's velocity and over time, backwards. The expected morphology 
of the remnant will therefore be of two expanding shells which will eventually join, producing yet 
another shock. At late times the shells will merge and become a single spherical shell. We made 
several runs, differing in the initial density, shape, initial velocity, angular width and the numerical 
coefficient in Eq. (1) of the ejecta as summarized in Table 1. In five runs the initial ejecta was 
shaped like a disk and in one case it was shaped like a sphere. As we show in Fig. 2, for ^ > 5 x 10^ 
the z positions of the shock and its maximal radius in the x — y plane, Vxy are indeed unaffected 
by these variations in initial conditions. The at which the shells touch depends on the exact 
initial conditions, but all other events (such as the ^ at which the remnant becomes spherical) are 
unaffected. Therefore in the following discussions we will consider only one of the runs. 

The mass inside the shell, defined as all mass with density below pq/2 also evolves almost 
linearly with ^l (Fig. 3) as 0.06(^o/c^)a*- In Fig- 4, we show the density contours as a function 
of n- The two shells touch and a shock forms along the equatorial plane between them at t ~ 
50 — 260(-E5i/pi)~^''^ yr when /x = 1 — 5 x 10^ depending on the initial conditions. The maximal z 
position of the shock and Txy can be fitted with a power law as shown in Fig. 5. The scaling relation 
for rxy is Txy = 0.07/x'^'^^(£J5i//3i)^/^ pc oc t^'^^ . The effective radius of the shell, {r'^yzY^^, scales as 

exactly the result expected for a Sedov-Taylor blast wave. The ratio Zmaxl^xy decreases in time 
staying always between 1 and 2 (see Fig 6) . Extrapolating we see that this ratio reaches a value of 
1 at ~ 9.6 X 10^. At this time the shock has a spherical shape with z = rxy ~ 12(£^5i/pi)^/^pc. 
Even then the shock will not be completely spherically symmetric as there would still be a ring of 
shocks around the "equator" where the shells have collided. 

During the whole run, total energy is conserved to within 1%. In Fig. 7 we show the total 
kinetic energy in the z direction, Exy the total kinetic energy in the x — y plane and Ei the total 
internal energy in the simulated volume. As the blobs interact with the ISM E^ and Exy increase 
at the expense of E^. Exy increases until fi ^ 5 x 10^ and then it remains constant at 0.22£'o. 
Ei increases by a factor of 200 to a value of 0.72^0) most of the increase occurs before /x ~ 2500. 
In the final configuration (^u ~ 2 x 10^) 72% of the energy is in internal energy, 22% is in kinetic 
energy in the x — y plane and only 6% remains in the kinetic energy in the z direction. 

Figures 8 and 9 depict the images of the remnant as a function of time and angles of inclination. 
We show images due to bremsstrahlung emission and synchrotron emission. The images are con- 
structed assuming that all the gas is optically thin in the relevant frequencies. The bremsstrahlung 
luminosity (Fig. 8) was calculated assuming that the volume cmissivity is proportional to p^e^/^ 
(Lang 1980). In calculating the synchrotron emissivity (Fig. 9) we assumed that both the magnetic 
field energy density and the energy of the relativistic electrons are proportional to the internal en- 
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ergy density of the gas with the proportionahty factors €b and €e respectively. We further assume 

that the relativistic electron number density is a power law in energy. Under these assumptions the 
volume emissivity is proportional to p^e^ (e.g. Shu 1991). In the late images there are two bright 
circles at the lines where the colliding blobs forming a hot shocked region. 

4. Discussion 

We follow the hydrodynamic evolution of two blobs colliding with the ISM. This scenario is a 
model for the intermediate time behavior of matter ejected from a central engine producing a GRB. 
Our results can be rescaled to fit a variety of initial energies and ISM densities. The two shells 
touch and a shock forms along the equatorial plane between them at f ~ 50 — 160(£^5i//9i)~^/^ yr 
when fi = 1 — 5 X 10^. We show that the late time remnant is insensitive to the exact initial 
morphology, angular width and density of the cjccta. Although initially the remnant may be highly 
non-spherical, the ratio between its height and radius will approach unity and it will eventually 
become spherical in shape after a time of ^ 5 x IO^eK^ p^^^^ yr when n>9.6x 10^. After this time 
it will be difficult to distinguish a GRB remnant from a SNR using the morphology. The expected 
number of non-spherical GRB remnants is, therefore, 5 x W~'^ ff^^ eI^^^ p^ per galaxy (using the 
observed present GRB rate (Schmidt 1999) of 10~^ yr~^ gal~^). 

The results of Gaensler (1998) show a tendency for the bilateral axis of the non-spherical SNRs 
to be aligned with the galactic plane. This presents strong evidence in favor of an extrinsic model 
for the origin of non-spherical SNRs. Using our results, we can propose an alternative explanation 
to the origin of some of the highly non-spherical SNRs and Hi shells. Instead of assuming an 
extrinsic model, namely spherical energy deposition into a non isotropic medium we propose an 
intrinsic model: non spherical energy deposition into an isotropic medium. Our GRB remnant 
model can explain non-spherical SNRs with a two shell morphology provided that p < 9.6 x 10^. 
One extreme and well studied case, the non-spherical SNR DEM L316 (Williams et al. 1997), has 
a distinct double shell morphology. The external model clearly fails here. However the ratio /i 
measured for DEM L316 is > 7.1 x 10^ and a a GRB remnant would already be spherical at this 
stage, suggesting that DEM L 316 is not a GRB remnant. This conclusion could be revised if for 
some reason the initial blobs stay confined for a much longer period (so that eq. 1 is not satisfied) 
or if the solution becomes radiative before becoming spherical and the radiative cooling slow down 
the evolution. 

We thank the anonymous referee for his helpful comments. This research was supported by a 
ISRAEL-US BSF grant. 
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run shape angular width density [pi] 



1 disc 10/3 6/5 

2 disc 1 4/3 

3 disc 2 3 

4 disc 1 3 

5 disc 7/12 3 

6 sphere 1/2 2 



Table 1: Initial parameters for the different runs. Initial density is in units of pi. 
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Fig. 1. — Time as a function of /x. The linear relation between time and /x is t ~ 
0.046/x(£^5i//9i)~^/^ yr. In this and all subsequent figures the results are presented for an initial 
energy of 10^^ erg and ISM density of 10~^^ gmcm~^. 
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Fig. 2. — Comparison of the various runs showing their similarities, (a) The z and r^y positions of 
the shocks (the z positions are the higher hnes). (b) The ratio zjvxy for ah the runs. 
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Fig. 3. — Mass inside the shell as a function of fi. 
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Fig. 4. — Density contours. The length scale is in pc. Contours are equal spaced at 
1.5po, 2po, 3.5po- The images are at ^ = 2.1 x 10^ 4.2 x 10^ 8.5 x 10^, 1.5 x 10^, 2.5 x 
10^, 4.5 X 10^, 7.6 X 10^, 1.3 x 10^^,2.2 x 10^ (left to right, top to bottom) 
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10^ 10^^ 10^ 

Fig. 5. — log-log plot of the z position of the blobs (crosses) and the maximum radius in the x — y 
plain of the shock rxy (circles) as a functions of /i. The solid lines are the best fit power laws for 
H > 2200 which are z oc /jP''^^ and Vxy oc /x*^'^^. The lines are extrapolated to the point where they 
cross at /X ~ 9.6 x 10^ 
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Fig. 7. — plot of total kinetic energy in the z direction (dashed line), total kinetic energy in the 
X — y direction (solid line) and total internal energy (dash-dotted line) . The energies are scaled by 
the initial kinetic energy Eq. 
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Fig. 8. — Images of the remnant, bremsstrahlung emission. The number above each image is the 
angle of incUnation in degrees. The images are shown at the same /x as the last 8 panels of figure 4. 
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Fig. 9. — Images of the remnant, synchrotron emission. The ji are the same as the last 8 panels in 
figure 8 



